###################
## Analysis
###################
library(tidyverse)
library(lme4)
library(ggeffects)
library(stargazer)
library(ggpubr)

# Set your working directory
setwd("~/Dropbox/2001_replication/draft/data")

# Load data
dat <- read.csv("output/rueda_rep.csv")

# Make sure that year and occupation are factor variables
dat$year <- as.factor(dat$year)
dat$euroesec <- as.factor(dat$euroesec)

################
## Table 1 
################

# Subgroup analysis

# Southern European Countries
fit_southern <- glmer(redpref1 ~ incomedif*forpop + age + gender +
                       education + religion + union + euroesec +
                       socgdp + year + (1 | cntry),
                     data = dat %>% filter(brncntr == 1, welfare == "southern"),
                     family = binomial(link = logit), nAGQ = 10)
summary(fit_southern)

# Non-Southern European Countries
fit_nonsouthern <- glmer(redpref1 ~ incomedif*forpop + age + gender +
                        education + religion + union + euroesec +
                        socgdp + year + (1 | cntry),
                      data = dat %>% filter(brncntr == 1, welfare != "southern"),
                      family = binomial(link = logit), nAGQ = 10)
summary(fit_nonsouthern)

# Nordic countries (social democratic)
fit_socdem <- glmer(redpref1 ~ incomedif*forpop + age + gender +
                       education + religion + union + euroesec +
                       socgdp + year + (1 | cntry),
                     data = dat %>% filter(brncntr == 1, welfare == "social_dem"),
                     family = binomial(link = logit), nAGQ = 10)
summary(fit_socdem)

# Continental European countries (corporatist)
fit_corp <- glmer(redpref1 ~ incomedif*forpop + age + gender +
                      education + religion + union + euroesec +
                      socgdp + year + (1 | cntry),
                    data = dat %>% filter(brncntr == 1, welfare %in% c("corporatist")),
                    family = binomial(link = logit), nAGQ = 10)
summary(fit_corp)

# Liberal countries
fit_lib <- glmer(redpref1 ~ incomedif*forpop + age + gender +
                    education + religion + union + euroesec +
                    socgdp + year + (1 | cntry),
                  data = dat %>% filter(brncntr == 1, welfare %in% c("liberal")),
                  family = binomial(link = logit), nAGQ = 10)
summary(fit_lib)

## Table 1
stargazer(fit_southern, fit_nonsouthern, fit_socdem, fit_corp, fit_lib, type = "latex")

###################
## Figures 3 & 4
###################

# using MCWR (to make Figure 4)
fit_d_mcwr <- glmer(redpref1 ~ incomedif*forpop*d_mcwr + age + gender + 
                         education + religion + union + euroesec +
                         socgdp + year + (1|cntry),
                       data = dat %>% filter(brncntr==1),
                       family = binomial(link = logit))
summary(fit_d_mcwr)

# For Robustness Check
fit_d_mcwr_fe <- glm(redpref1 ~ incomedif*forpop*d_mcwr + age + gender + 
                      education + religion + union + euroesec +
                      socgdp + year + cntry,
                    data = dat %>% filter(brncntr==1),
                    family = binomial(link = logit))

# Function for plotting predicted probabilities
pred_plot <- function(pred_dat, ylab = "Predicted Probabilities of Support for Redistribution"){
  if ("facet" %in% colnames(pred_dat)){
    p <- pred_dat %>% 
      mutate(group = as.numeric(as.character(group))) %>%
      mutate(Group = ifelse(group==min(group),
                            "Low Foreign-born Population",
                            "High Foreign-born Population")) %>%
      ggplot(aes(x=x, y=predicted, group = Group)) +
      geom_line(aes(col=Group), size = 0.7) +
      geom_ribbon(aes(ymin=conf.low, ymax=conf.high, fill = Group), alpha = 0.2) +
      theme_bw() + theme(legend.position = "bottom") +
      scale_color_manual(values=c("orangered3", "steelblue4")) +
      xlab("Distance to Mean Income (10,000s PPP 2010 dollars)") + ylab(ylab) +
      facet_wrap(~facet)
    return(p)
  } else {
  p <- pred_dat %>% 
    mutate(group = as.numeric(as.character(group))) %>%
    mutate(Group = ifelse(group==min(group),
                          "Low Foreign-born Population",
                          "High Foreign-born Population")) %>%
    ggplot(aes(x=x, y=predicted, group = Group)) +
    geom_line(aes(col=Group), size = 0.7) +
    geom_ribbon(aes(ymin=conf.low, ymax=conf.high, fill = Group), alpha = 0.2) +
    theme_bw() + theme(legend.position = c(0.25, 0.2)) +
    scale_color_manual(values=c("orangered3", "steelblue4")) +
    xlab("Distance to Mean Income (10,000s PPP 2010 dollars)") + ylab(ylab)
  return(p)
  }
}

###############
# Figure 3
###############
pred_southern <- ggpredict(fit_southern, c("incomedif [-2.5:9.2]","forpop[5.3, 14.7]"))
pred_plot(pred_southern)

pred_socdem <- ggpredict(fit_socdem, c("incomedif [-2.5:9.2]","forpop[5.3, 14.7]"))
pred_plot(pred_socdem)

pred_corp <- ggpredict(fit_corp, c("incomedif [-2.5:9.2]","forpop[5.3, 14.7]"))
pred_plot(pred_corp)

pred_lib <- ggpredict(fit_lib, c("incomedif [-2.5:9.2]","forpop[5.3, 14.7]"))
pred_plot(pred_lib)

pred_southern$facet <- "Southern"
pred_socdem$facet <- "Nordic"
pred_corp$facet <- "Continental"
pred_lib$facet <- "Liberal"
pred_regime <- rbind(pred_southern, pred_socdem, pred_corp, pred_lib)

# Figure 3
pred_plot(pred_regime)

###############
# Figure 4
###############
pred_mcwr <- ggpredict(fit_d_mcwr, c("incomedif [-2.5:9.2]","forpop[5.3, 14.7]", "d_mcwr"))
pred_mcwr$facet <- ifelse(pred_mcwr$facet==0, "Low Migrant Care Worker Ratio", "High Migrant Care Worker Ratio")

# Figure 4
pred_plot(pred_mcwr)

#####################################
## [Figure 5 & 6] Micro-level mechanisms
#####################################

# Child care services for working parents, governments' responsibility
fit_childcare <- lmer(childcare ~ incomedif*forpop*d_mcwr + age + gender + education + religion + union + euroesec +
                        socgdp + (1 | cntry),
                      data = dat %>% filter(brncntr == 1))
summary(fit_childcare)

# Standard of living for the unemployed, governments' responsibility 
fit_unemp <- lmer(unemployed ~ incomedif*forpop*d_mcwr + age + gender + education + religion + union + euroesec +
                      socgdp + (1 | cntry), data = dat %>% filter(brncntr == 1))
summary(fit_unemp)

unique(dat$cntry[dat$d_mcwr==0 & dat$essround==4])

# Figure 5
pred_childcare <- ggpredict(fit_childcare, c("incomedif [-2.5:9.2]","forpop[5.3, 14.7]", "d_mcwr")) 
pred_childcare$facet <- ifelse(pred_childcare$facet==0, "Low Migrant Care Worker Ratio", "High Migrant Care Worker Ratio")
pred_plot(pred_childcare)

# Figure 6
pred_unemp <- ggpredict(fit_unemp, c("incomedif [-2.5:9.2]","forpop[5.3, 14.7]", "d_mcwr"))
pred_unemp$facet <- ifelse(pred_unemp$facet==0, "Low Migrant Care Worker Ratio", "High Migrant Care Worker Ratio")
pred_plot(pred_unemp)

#####################
# Appendix
#####################
# Interaction models
fit_welfare <- glmer(redpref1 ~ incomedif*forpop*welfare + age + gender +
                       education + religion + union + euroesec +
                       socgdp + year + (1 | cntry),
                     data = dat %>% filter(brncntr == 1),
                     family = binomial(link = logit), nAGQ = 10)
summary(fit_welfare)

# Table A1
stargazer(fit_welfare, fit_d_mcwr, fit_d_mcwr_fe, type = "text")

# Figure: Interaction Model
pred_welfare <- ggpredict(fit_welfare, c("incomedif [-2.5:9.2]","forpop[5.3, 14.7]","welfare"))
pred_welfare$facet <- recode(pred_welfare$facet, "corporatist" = "Continental", "liberal" = "Liberal",
                             "social_dem" = "Nordic", "southern" = "Southern")
pred_plot(pred_welfare)

# Allow immigrants of different race/ethnic group from majority
fit_imm1 <- lmer(imm_etn ~ incomedif*forpop*d_mcwr + age + gender + education + religion + union + euroesec +
                   socgdp + (1 | cntry), data = dat %>% filter(brncntr == 1))

# Allow immigrants from poorer countries outside Europe
fit_imm2 <- lmer(imm_pcntr ~ incomedif*forpop*d_mcwr + age + gender + education + religion + union + euroesec +
                   socgdp + (1 | cntry), data = dat %>% filter(brncntr == 1))

# Figure 7
pred_imm1 <- ggpredict(fit_imm1, c("incomedif [-2.5:9.2]","forpop[5.3, 14.7]", "d_mcwr")) 
pred_imm1$facet <- ifelse(pred_imm1$facet==0, "Low Migrant Care Worker Ratio", "High Migrant Care Worker Ratio")
pred_plot(pred_imm1, ylab = "Support for Immigration")

# Figure 8
pred_imm2 <- ggpredict(fit_imm2, c("incomedif [-2.5:9.2]", "forpop[5.3, 14.7]", "d_mcwr")) 
pred_imm2$facet <- ifelse(pred_imm2$facet==0, "Low Migrant Care Worker Ratio", "High Migrant Care Worker Ratio")
pred_plot(pred_imm2, ylab = "Support for Immigration")

# Save the estimation results
save(fit_southern, fit_nonsouthern, fit_socdem, fit_corp, fit_lib, # Table 1
     fit_welfare, fit_d_mcwr, # Appendix Table 2
     fit_childcare, fit_unemp, 
     fit_imm1, fit_imm2,
     pred_plot,
     file = "output/final_results.RData") 
